Learn R Programming

frailtypack (version 2.3)

frailtyPenal for Joint frailty models: Fit Joint Frailty model for recurrent and terminal events using semiparametric penalized likelihood estimation or a parametric estimation

Description

Fit a joint either with gamma or log-normal frailty model for recurrent and terminal events using a penalized likelihood estimation on the hazard function or a parametric estimation. Right-censored data is allowed. Left-truncated data and stratified analysis are not possible. Joint frailty models allow studying, jointly, survival processes of recurrent and terminal events, by considering the terminal event as an informative censoring. There is two kinds of joint frailty models that can be fitted with frailtyPenal : - The first one (Rondeau et al. 2007) includes a common frailty term to the individuals $(\omega_i)$ for the two rates which will take into account the heterogeneity in the data, associated with unobserved covariates. The frailty term acts differently for the two rates ( $\omega_i$ for the recurrent rate and $\omega_i^{\alpha}$ for the death rate). The covariates could be different for the recurrent rate and death rate. For the $j^{th}$ recurrence $(j=1,...,n_i)$ and the $i^{th}$ subject $(i=1,...,G)$, the joint gamma frailty model for recurrent event hazard function $r_{ij}(.)$ and death rate $\lambda_i(.)$ is : $$\left{ \begin{array}{ll} r_{ij}(t|\omega_i)=\omega_ir_0(t)\exp(\bold{\beta_1^{'}Z_i(t)}) & \mbox{(Recurrent)} \ \lambda_i(t|\omega_i)=\omega_i^{\alpha}\lambda_0(t)\exp(\bold{\beta_2^{'}Z_i(t)}) & \mbox{(Death)} \ \end{array} \right.$$ where $r_0(t)$ (resp. $\lambda_0(t)$) is the recurrent (resp. terminal) event baseline hazard function, $\bold{\beta_1}$ (resp. $\bold{\beta_2}$) the regression coefficient vector, $\bold{Z_i(t)}$ the covariate vector. The random effects of frailties $\omega_i\sim\bold{\Gamma}(\frac{1}{\theta},\frac{1}{\theta})$ and are iid. The joint log-normal frailty model will be : $$\left{ \begin{array}{ll} r_{ij}(t|\eta_i)=r_0(t)\exp(\eta_i+\bold{\beta_1^{'}Z_i(t)}) & \mbox{(Recurrent)} \ \lambda_i(t|\eta_i)=\lambda_0(t)\exp(\alpha \eta_i+\bold{\beta_2^{'}Z_i(t)}) & \mbox{(Death)} \ \end{array} \right.$$ where $$\eta_i\sim N(0,\sigma^2)$$ - The second one (Rondeau et al. 2011) is quite similar but the frailty term is common to the individuals from a same group. This model is useful for joint modelling two clustered survival outcomes. In this case, j is for the subject and i for the cluster. $$\left{ \begin{array}{ll} r_{ij}(t|u_i)=u_ir_0(t)\exp(\bold{\beta_1^{'}Z_{ij}(t)}) & \mbox{(Time to event)} \ \lambda_{ij}(t|u_i)=u_i^{\alpha}\lambda_0(t)\exp(\bold{\beta_2^{'}Z_{ij}(t)}) & \mbox{(Death)} \ \end{array} \right.$$ In case of a log-normal distribution of the frailties, we will have : $$\left{ \begin{array}{ll} r_{ij}(t|v_i)=r_0(t)\exp(v_i+\bold{\beta_1^{'}Z_{ij}(t)}) & \mbox{(Time to event)} \ \lambda_{ij}(t|v_i)=\lambda_0(t)\exp(\alpha v_i+\bold{\beta_2^{'}Z_{ij}(t)}) & \mbox{(Death)} \ \end{array} \right.$$ where $$v_i\sim N(0,\sigma^2)$$ This joint frailty model can also be applied to clustered recurrent events and a terminal event (example on "readmission" data below).

Arguments

formula
a formula object, with the response on the left of a $\texttildelow$ operator, and the terms on the right. The response must be a survival object as returned by the 'Surv' function like in survival package. To fit a joint
formula.terminalEvent
a formula object, only requires terms on the right to indicate which variables are modelling the terminal event.
data
a 'data.frame' in which to interpret the variables named in the 'formula' and 'formula.terminalEvent'.
Frailty
Logical value. Is model with frailties fitted? If so, variance of frailty parameter is estimated. The default is FALSE
joint
Logical value. Is joint model fitted? If so, 'formula.terminalEvent' is required. The default is FALSE
recurrentAG
Logical value. Is Andersen-Gill model fitted? If so indicates that recurrent event times with the counting process approach of Andersen and Gill is used. This formulation can be used for dealing with time-dependent cov
cross.validation
Logical value. Is cross validation procedure used for estimating smoothing parameter in the penalized likelihood estimation? If so a search of the smoothing parameter using cross validation is done, with kappa1 as the seed.
n.knots
integer giving the number of knots to use. Value required in the penalized likelihood estimation. It corresponds to the (n.knots+2) splines functions for the approximation of the hazard or the survival functions. Number o
kappa1
positive smoothing parameter in the penalized likelihood estimation. The coefficient kappa of the integral of the squared second derivative of hazard function in the fit (penalized log likelihood). To obtain a good value for
kappa2
positive smoothing parameter in the penalized likelihood estimation for the terminal event rate. To obtain a good value for kappa2, a solution is to fit the corresponding Cox model using cross-validation (See cross.validation)
maxit
maximum number of iterations for the Marquardt algorithm. Default is 350
hazard
Type of hazard functions: "Splines" for semiparametric hazard functions with the penalized likelihood estimation, "Piecewise-per" for piecewise constant hazard function using percentile, "Piecewise-equi" for piecewise constant
nb.int1
Number of intervals (between 1 and 20) for the recurrent parametric hazard functions ("Piecewise-per", "Piecewise-equi").
nb.int2
Number of intervals (between 1 and 20) for the death parametric hazard functions ("Piecewise-per", "Piecewise-equi").
RandDist
Type of random effect distribution: "Gamma" for a gamma distribution, "LogN" for a log-normal distribution. Default is "Gamma".

Value

  • Parameters estimates of a joint frailty model, more generally a 'fraityPenal' object. Methods defined for 'frailtyPenal' objects are provided for print, plot and summary. The following components are included in a 'frailtyPenal' object for Joint frailty models.
  • bsequence of the corresponding estimation of the splines coefficients, the random effects variances, the power of the frailties and the regression coefficients.
  • alphathe coefficient associated with the frailty parameter in the terminal hazard function.
  • callThe code used for fitting the model.
  • coefthe regression coefficients.
  • cross.ValLogical value. Is cross validation procedure used for estimating the smoothing parameters in the penalized likelihood estimation?
  • formulathe formula part of the code used for the model.
  • groupsthe maximum number of groups used in the fit.
  • kappaA vector with the smoothing parameters in the penalized likelihood estimation corresponding to each baseline function as components.
  • lammatrix of hazard estimates and confidence bands.
  • lam2the same value as lam for the second stratum.
  • loglikPenalthe complete marginal penalized log-likelihood in the semiparametric case.
  • loglikthe marginal log-likelihood in the parametric case.
  • nthe number of observations used in the fit.
  • n.eventsthe number of events observed in the fit.
  • n.knotsnumber of knots for estimating the baseline functions.
  • n.iternumber of iterations needed to converge.
  • survmatrix of baseline survival estimates and confidence bands.
  • surv2the same value as surv for the the second stratum.
  • thetavariance of the gamma frailty parameter $(\bold{Var}(\omega_i))$ or $(\bold{Var}(u_i))$
  • sigma2variance of the log-normal frailty parameter $(\bold{Var}(\eta_i))$ or $(\bold{Var}(v_i))$
  • varHthe variance matrix of all parameters before positivity constraint transformation (theta, the regression coefficients and the spline coefficients). Thenafter, the delta method is needed to obtain the estimated variance parameters.
  • varHIHthe robust estimation of the variance matrix of all parameters (theta, the regression coefficients and the spline coefficients).
  • x1vector of times where both survival and hazard functions for the recurrent events are estimated. By default seq(0,max(time),length=99), where time is the vector of survival times.
  • x2vector of times for the second stratum (see x1 value).
  • type.of.hazardType of hazard functions (0:"Splines", "1:Piecewise", "2:Weibull").
  • type.of.PiecewiseType of Piecewise hazard functions (1:"percentile", 0:"equidistant").
  • nbintervRNumber of intervals (between 1 and 20) for the recurrent parametric hazard functions ("Piecewise-per", "Piecewise-equi").
  • nbintervDCNumber of intervals (between 1 and 20) for the death parametric hazard functions ("Piecewise-per", "Piecewise-equi").
  • nparnumber of parameters.
  • nvarA vector with the number of covariates of each type of hazard function as components.
  • nvarRecnumber of recurrent explanatory vairables.
  • nvarEndnumber of death explanatory vairables.
  • noVar1indicator of recurrent explanatory vairables.
  • noVar2indicator of death explanatory vairables.
  • LCVthe approximated likelihood cross-validation criterion in the semiparametric case (with H minus the converged hessien matrix, and l(.) the full log-likelihood.$$LCV=\frac{1}{n}(trace(H^{-1}_{pl}H) - l(.))$$
  • AICthe Akaike information Criterion for the parametric case.$$AIC=\frac{1}{n}(np - l(.))$$
  • n.knots.tempinitial value for the number of knots.
  • shape.weibshape parameters for the weibull hazard function.
  • scale.weibscale parameters for the weibull hazard function.
  • martingale.resmartingale residuals for each cluster (recurrent).
  • martingaledeath.resmartingale residuals for each cluster (death).
  • frailty.predempirical Bayes prediction of the frailty term.
  • linear.predlinear predictor: uses "Beta'X + log w_i" in the joint gamma frailty models, otherwise uses "Beta'X + eta_i" in log-normal
  • lineardeath.predlinear predictor for the terminal part form the joint model: "Beta'X + alpha.log w_i" for gamma, "Beta'X + alpha.eta_i" .
  • global_chisqa vector with the values of each multivariate Wald test.
  • dof_chisqa vector with the degree of freedom for each multivariate Wald test.
  • global_chisq.testa binary variable equals to 0 when no multivariate Wald is given, 1 otherwise.
  • p.global_chisqa vector with the p_values for each global multivariate Wald test.
  • names.factorNames of the "as.factor" variables.
  • Xlevelsvector of the values that factor might have taken for the recurrent part.
  • contraststype of contrast for factor variable for the recurrent part.
  • Xlevels2vector of the values that factor might have taken for the death part.
  • contrasts2type of contrast for factor variable for the death part.

synopsis

frailtyPenal(formula, formula.terminalEvent, data, Frailty = FALSE, joint = FALSE, recurrentAG = FALSE, cross.validation = FALSE, n.knots, kappa1, kappa2, maxit = 350, hazard = "Splines", nb.int1, nb.int2, RandDist = "Gamma")

Details

The estimated parameter are obtained by maximizing the semiparametric penalized loglikelihood estimation or a parametric estimation using the robust Marquardt algorithm (Marquardt, 1963) which is a combination between a Newton-Raphson algorithm and a steepest descent algorithm. When frailty parameter is small, numerical problems may arise. To solve this problem, an alternative formula of the penalized log-likelihood is used (see Rondeau, 2003 for further details). Cubic M-splines of order 4 are used for the hazard function, and I-splines (integrated M-splines) are used for the cumulative hazard function. The inverse of the hessian matrix is the variance estimator and to deal with thepositivity constraint of the variance component, a squared transformation is used and the Standard errors are computed by the $\Delta$-method (Knight & Xekalaki, 2000). The smooth parameter can be chosen by maximizing a likelihood cross validation criterion (Joly and other, 1998). The integrations in the full log likelihood were evaluated using Gaussian quadrature. Laguerre polynomials with 20 points were used to treat the integrations on $[0,\infty[$. INITIAL VALUES The splines and the regression coefficients are initialized to 0.5. The program fits an adjusted Cox model to have new initial values for the regression and the splines coefficients. The variance of the frailty term $\theta$ and the coefficient $\alpha$ associated in the death hazard function are initialized to 1. Then, it fits a joint frailty model.

References

V. Rondeau, Y. Mazroui and J. R. Gonzalez (2012). Frailtypack: An R package for the analysis of correlated survival data with frailty models using penalized likelihood estimation or parametric estimation. Journal of Statistical Software 47, 1-28. V. Rondeau, J.P. Pignon, S. Michiels (2011). A joint model for the dependance between clustered times to tumour progression and deaths: A meta-analysis of chemotherapy in head and neck cancer. Statistical methods in medical research 897, 1-19. V. Rondeau, S. Mathoulin-Pellissier, H. Jacqmin-Gadda, V. Brouste, P. Soubeyran (2007). Joint frailty models for recurring events and death using maximum penalized likelihood estimation:application on cancer events. Biostatistics 8,4, 708-721. V. Rondeau, D. Commenges, and P. Joly (2003). Maximum penalized likelihood estimation in a gamma-frailty model. Lifetime Data Analysis 9, 139-153. D. Marquardt (1963). An algorithm for least-squares estimation of nonlinear parameters. SIAM Journal of Applied Mathematics, 431-441.

See Also

summary.jointPenal, print.jointPenal, plot.jointPenal, readmission, terminal, cluster

Examples

Run this code
### Joint model (recurrent and terminal events) with 2 covariates ###

data(readmission)

### Joint Frailty model ###

## Gap-time ##
modJoint.gap <- frailtyPenal(Surv(time,event)~cluster(id)+sex+dukes+charlson+
                  terminal(death),formula.terminalEvent=~sex+dukes+charlson,
                  data=readmission,n.knots=14,kappa1=9.55e+9,kappa2=1.41e+12,
                  Frailty=TRUE,joint=TRUE,recurrentAG=FALSE)

## Calendar time ##
modJoint.calendar <- frailtyPenal(Surv(t.start,t.stop,event)~cluster(id)+sex
              +dukes+charlson+terminal(death),formula.terminalEvent=~sex
              +dukes+charlson,data=readmission,n.knots=10,kappa1=9.55e9,
              kappa2=1.41e12,Frailty=TRUE,joint=TRUE,recurrentAG=TRUE)

### Joint Frailty model with log-normal distribution ###

modJoint.log <- frailtyPenal(Surv(t.start,t.stop,event)~cluster(id)+sex
              +dukes+charlson+terminal(death),formula.terminalEvent=~sex
              +dukes+charlson,data=readmission,n.knots=10,kappa1=9.55e9,
              kappa2=1.41e12,Frailty=TRUE,joint=TRUE,recurrentAG=TRUE,
              RandDist="LogN")

### Joint Frailty model for clustered data ###

# here is a generated cluster (31 clusters of 13 subjects)
readmission <- transform(readmission,group=id%%31+1)

joi.clus <- frailtyPenal(Surv(t.start,t.stop,event)~cluster(group)+
              num.id(id)+dukes+charlson+sex+chemo+terminal(death),
              formula.terminalEvent=~dukes+charlson+sex+chemo,
              data=readmission,Frailty=TRUE,joint=TRUE,recurrentAG=TRUE,
              n.knots=10,kappa1=2.11e+08,kappa2=9.53e+11)

Run the code above in your browser using DataLab